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Abstract 

This paper covers the material of our two talks. We describe a series of projects based upon 
perturbative expansions to follow the gravitational evolution of the one point probability distribu- 
tion functions (PDFs) for the density contrast field and for the divergence of the corresponding 
velocity field. The Edgeworth expansion greatly simplifies the problem. Indeed, the PDF can be 
described through several low order moments, with the first non-trivial contribution coming from 
the skewness, or third moment. 

Gaussian initial conditions imply that the skewness of the smoothed density field is proportional 
to the square of the variance. The constant of proportionality is then a function of the logarithmic 
slope of the variance versus scale, and is otherwise independent of scale. This constant is also 
practically independent of the density parameter SI, and has very nearly the same value when 
distances are estimated by redshifts. To show the latter property, we also briefiy discuss our 
Lagrangian approach to perturbation theory. Finally, the observed skewness does depend, of course, 
on the relation between the distribution of light and mass; we describe the case of local biasing. 

For the divergence of the velocity field, the skewness of its PDF is also proportional to the 
square of its variance, but the ratio now strongly depends on the value of SI. This offers a new 
method to determine the density parameter which does not involve comparisons with the density 
field; it should be independent of assumptions concerning the way light traces mass. 

Comparisons with numerical simulations show that these theoretical results can be applied to 
the real non-linear universe. We thus proceed to compare the density and velocity field skewness 
with observational data (IRAS and POTENT). The results are compatible with gaussian initial 
conditions, and local biasing. 



1 The PDF as a Measure of Large Scale Structures 

Many approaches can be used to characterize statistically the observed Large Scale Structures in the 
Universe. Here we focus on the one-point distribution function, or PDF. Empirically, this approach 
was first applied by Hubble (1934). Let Q(r) be a physical quantity whose value depends on the 
location. Let Qi{r) be the smoothed value of Q on scale I, i.e., Qi{r) is obtained by convolving the 
scalar field Q{r) by a smoothing window Wi of characteristic scale I (Qi(r) = Q(s) * Wi(r — s) = 



/ d'^s Q(s) Wi(r — s)). The probability that Qi lies between Q and Q + dQ is P(Q,i) dQ, P being its 
PDF. 

If Q = iV is a number of galaxies, and is a top-hat window, then the PDF gives the probability 
of finding N galaxies in a sphere (in 3D) of size i. In the following, Q will be either the mass density 
contrast field S = p/]) — 1, or the divergence of the associated velocity field in units of the Hubble 
constant, 9 = V -v/H; we will consider the two most commonly used types of windows, a top-hat and 
a gaussian one. Later on, we shall discuss the relation of the mass distribution to the distribution of 
luminous galaxies. 

It is frequently assumed that, very early in the history of the Universe, the density contrast 
e = 8{t = ti) could be taken as normally distributed, i.e., 



P(e, oc exp 



where ct'^{1) = (e|) stands for the initial variance of the density contrast field smoothed on scale I. It 
is assumed to be small at all scales, Ui <^1. The second moment is the only quantity we need to know 
in order to fully characterize such a field. Indeed all odd moments of a gaussian PDF are zero, while 
the 2ra-th moments are proportional to the n power of the variance, (e^") oc af'^ (the corresponding 
specification in Fourier space is given in (7)). In the following, we shall refer to such a state as gaussian 
initial conditions. 

The P(S,i) moments are simply related to other statistical quantities. For instance, 
(ST) = {6(ri)*Wi(r-ri)6(r2)*Wi(r-r2) ... 6(r^)*Wi(r-r^)) 
d^ri d^r2 . . . d^r^ Wi(ri)Wi(r2) . . .Wi(r^) F„(ri, . . . r„). 



that is the m-th order moment is the volume average (with weights Wi) of the full m-point correlation 
function, Fmi'ri,r2 . . .r^) = (^(ri)^(r2) . . .S(rm)) (since the overall homogeneity and isotropy of the 
system implies = (^(ri — r)^(r2 — r) . . .^(r — Fm)))- These correlation functions have been the 
most widely used indicators so far to describe the small-scale galaxy distribution; they give the proba- 
bility of finding m galaxies in infinitesimal volumes around ri, r2 . . . r^; only the lowest orders (m < 4) 
are well-known (and the reduced parts of F2, F3, and F4 are usually denoted (respectively) by = 

C = F3, and ?7(xi,X2,X3,X4) = F4 - F2(xi, X2)-F2(X3, X4) - F2(xi, X3)F2(X2, X4) - F2(xi, X4)F2(X2, X3); 

see LSS for further details). Also, if the m-th order moment behaves as a power law of scale i, it can 
be used to obtain the Renyi indices of multifractal analysis ((^™) oc see e.g., [2] for 

further details). The (one-point) PDF gives only a partial statistical description of a system. A com- 
plete description calls for an infinite hierarchy of m-point PDFs (the m-point PDF is the probability 
of finding ^^(1), Si(2),... Si(m) at points ri, r2, . . ., fm)- Still we shall see that this simple object 
contains a wealth of useful information. 

Our first task will be to derive the properties of the PDF today, under the assumption of a time- 
evolution under the sole infiuence of the gravitational instability acting on gaussian initial conditions. 
Under the infiuence of gravity, underdense regions become even more underdense (although not indef- 
initely, since the density is bounded from below), while positive density enhancements tend to grow 
without bound. Clearly the symmetry of the distribution cannot be maintained, and the PDF becomes 
skewed, i.e., (^■^) departs from zero. The distribution also develops a non-zero kurtosis (^'*) — 3 (^^) . 
In the limit of a small variance, (^^) ^ 1, the development of skewness and kurtosis are the most 
important effects. 

In the following, we first recall how one obtains the lowest order moments of the PDF in the weakly 
non-linear regime in Eulerian perturbation theory. We then address the question of the i7 dependence 
of the results, and of smoothing. We show next the relationship between the low order moments and 
the overall shape of the PDF, in the case of gaussian initial conditions. We proceed by discussing 
the overall validity range of this approach. We also show that Lagrangian perturbation theory is well 
suited to compute the effect on the density field of using redshifts as distance estimates. We conclude 
by a discussion of the implications of these results by comparing with observational data. 



2 Perturbation Theory 



Here we follow the standard Eulerian approach, as described in Peebles' book [36, hereafter LSS] in 
§18. It is assumed that matter may be described as a non-relativistic pressureless fluid embedded in 
a Friedman-Lemaitre model with zero cosmological constant. 



2.1 Method 



One starts from the standard equations for a perfect fluid, i.e., the continuity equation, Euler equation, 
and Poisson equation. We write the perturbative expansion for the density contrast field as 

where x are the comoving coordinates, and t is the cosmological time. A similar expansion is performed 
on the proper peculiar velocity v. The term is the linear order solution of the hydrodynamic 
equations of motion in comoving coordinates, 

^^^^ = _D(i)e(x) + decaying mode , (1) 

where e depends on the spatial coordinates alone and D(t) is the standard growing mode (LSS, 
eq. [11.16]). The term = 0(^(^))^ = O(e^) is the solution of the equations of motion with quadratic 
nonlinearities included iteratively by using ^^^^ as source terms (as in LSS, §18). The dominant mode 
in the second order solution is 



|(1 + k)s^ + Ve • + (^ - 



where V = d/dx while <?!>(x) = — / (i"^a;'e(x')/47r|x — x'| is the linearized Newtonian gravitational 
potential, and 

3 

= XI ^i*^ Tap = (ISapV^ -VaVp)(t). 

a,fj=l 

Here a, (3 are the indices of the spatial coordinate components, x = Apart from a multiplicative 

factor, Tafj is the linear order peculiar velocity shear tensor (LSS, eq. [14.12]). The parameter K(t) 
is a slowly varying function of cosmological time. Bouchet et al. (1993, Ref. [8]) show that it is well 
approximated by 

^ l,^-'"' , (3) 

in the range 0.05 > i7 > 3 (the accuracy of this approximation is then better than 0.4%). For i7 = 0, 
K = |. The exact expression for «(!)), valid in the entire range i7 > 0, is given in [8]). For i7 = 1, 
K = and we recover the well known Einstein-de Sitter solution {e.g., LSS, eq. [18.8] or Ref. [25]). 

2.2 Skewness for the Density Field 



Now we can calculate the gravitationally induced skewness under the assumption that e is a random 
gaussian field. The lowest order terms in the series for {6'^) are 

The linear solution (1) implies that the first term is D(t)'^ times the initial skewness, which is zero for 
gaussian initial conditions. The second order solution (2) shows that the second term involves (e'*), 
which is oc a"^ for gaussian initial conditions. Thus the skewness ratio 



(5) 



is a constant versus scale, and one finds (ref. [8]) up to terms of order (since (e^) = 0), 



Ss = 4 + 4k(1)) ~ y + ^(11-2/63 _ 1) . (6) 

The ~ sign above applies to the range of applicability of equation (3). The first term of this equation, 
34/7, had been obtained by Peebles more than a decade ago (LSS, §18). The weak i7-dependence of 
the full expression shows that nearly all the i7-dependence of the skewness (S'^) comes from that of 
the square of the variance. It simply refiects the fact that the second order growth rate is nearly equal 
to D(ty, as can be seen from (2). 

2.3 Smoothing 

To make contact with observables, we want the skewness of the density field Si, when smoothed with 
either a top-hat or a gaussian (spherically symmetric) window, which satisfies 

„2 j3„ 



W{x)d^x = l, and j W{x) x^ d^x = r . 

The first equation insures a proper normalization to unity, while the second requires the effective 
half-width i to be finite. The top hat case is appropriate for comparisons with the observed frequency 
distribution of galaxies. (In that case, discreteness corrections must be taken into account before 
comparing with the theory, e.g.^ (^(N/N — 1)^^ = 1/N + (^^) if one adopts the Poisson model, see 
e.g., LSS §33. A gaussian window, on the other hand, removes these discreteness fiuctuations, since 
it does not have sharp edges.) Here we recall the main results of [27]. 

The calculations are most conveniently done in Fourier space. For the initial field e, the Fourier 
components are given by 

ek = (^2Tr f/2 j ^^^"^ exp(ik • x) d^x . 
For gaussian initial conditions, one has 

(ekEk,) = ^£,(k V)P{k), (ekek'^q) = 0, 

(7) 

(ekEk'eqeq') = P{k)P{q)f>D{<\ + q')^B(k + k') cycl. (two terms) , 

where So is the Dirac delta, and P{k) is the (initial) power spectrum. Now we can use the above 
expressions for the first few moments, together with equations (1), (2), and (4) to derive 6*3 to lowest 
non- vanishing order, 

^3 = / 00 P{k) P(k') Wk Wy PF|k_k'| r(k, k') + 0{a^) . (8) 

Here r(k, k') stands for r(k, k') = 4 4k(1)) - 6i^{k/k') + [2- 4k(1))] P2(m), where = k • k'/kk', 
and P2 is a Legendre polynomial. If there is no smoothing (I^k = 1)? the dipole and quadrupole 
terms integrate to zero, and one simply recovers (6). On the other hand, as soon as one introduces 
smoothing, the result does depend on the initial power spectrum P{k). 

Let us assume that P{k) oc A;". Then, for a top-hat smoothing, and —3 < < 1, the equation (8) 
yields after painful calculations the simple result 

53 = 4 + 4K(l))-(3 + n), (9) 

A similar feat can also be accomplished for a gaussian window, although the result is slightly less 
simple (see figure 4). Actually, a careful inspection of the expression (8) shows that 5*3 should only 
depends on the effective (logarithmic) slope of the power spectrum at the smoothing scale. This was 
confirmed in [27] by comparing numerical integration for a CDM power spectrum with a prediction 



using (9). Our result was recently generalized by Bernardeau in [5], who showed that for an arbitrary 
gaussian field, smoothed with a top hat filter on a scale I, 



9 log (^2 \ 

53(^) = 4 + 4k(1))-7,, with 7, = -^Ai_2. (10) 
For a pure power law Pik), we have 7 = 3 + ra, in agreement with (9). 
2.4 Skewness for the Velocity Field 

The continuity equation gives the divergence of the velocity field in terms of the time derivative of the 
density contrast. Fet us call T3 the skewness ratio for the divergence of the velocity field in units of 
the Hubble constant (9 = V ■ v/H). One then finds (ref. [6]; see also equation (28) in §2.7.2), for a 
top-hat smoothing. 



26 

y - 7^ 



The asymmetry in the distribution of 9 is directly related to the asymmetry in S: voids and clusters 
in the mass distribution correspond to sources (9 > 0) and sinks (9 < 0) in the velocity field. 

2.5 Shape of the PDFs and Edgeworth Expansion 

We now wish to examine how gravitational instability drives a PDF away from its initial state, which 
we assume to be gaussian. We start by the density field PDF, P(S,i). We need to introduce the Gram- 
Charlier expansion, which allows one to reconstruct the PDF from its moments. Then we rearrange 
the Gram-Charlier series by collecting all terms of the same order. The result is the proper asymptotic 
expansion of the PDF in powers of ai, that we introduced in [28]. 

Fet us consider p(//), the PDF of the density field in terms of the standardized random variable 
z/ = Sf/ai. Fet us also introduce <?!>(/^) = (27r)~^/^ exp( — z/^/2), a gaussian (or Normal) PDF. Since we 
want to describe an evolution from gaussian initial conditions, it makes sense to consider an expansion 
of p(z/) in terms of <?!>(/^) and its derivatives. The Gram-Charlier series (Cramer 1946 and references 
therein) provides such an expansion: 

K//) = CO H'^) + ^ <^(^ V) + 1 ^^'H'^) + • • • , (12) 
where are constant coefficients. Superscripts denote derivatives with respect to z/: 

^'•'^\^)^^ = {-irHU^)cp{v\ (13) 

where Hm is the Hermite polynomial of degree m. 

The Hermite polynomials satisfy orthogonality relations (e.g. Abramowitz & Stegun 1964): 



H^{v) Hp{v) (]){v) dv 



0, if m 7^ p ; 
m ! otherwise. 



Therefore, multiplying both sides of equation (12) by Hm and integrating term by term yields 

/•oo 

cm = {-ir / H^{v)p{v)dv. (14) 



Equation (14) gives cq = 1, ci = C2 = 0, while for the next four coefficients in the series we obtain 

c„ = (-1)" S.mcjf-'^ , for 3 < m < 5 ; cq = Seaj + 105|(t| . (15) 

Thus the Sm{£) have both a dynamic and a static application: they describe the time evolution 
of moments of the PDF at a fixed smoothing scale I, and they also describe the relation between 



moments of the PDF at a fixed time on different smoothing scales. In the latter case, one must also 
include the scale-dependence of the Sm if the initial power spectrum is not scale-free. 

We have just seen that perturbation theory and gaussian initial conditions imply that (^|) oc aj, 
and 5*3 is therefore an "order unity" quantity when ai is the "small" parameter. The same is true 
for all remaining reduced moments, Sm = 0(1) for all m (Bernardeau 1992, Fry 1994). For our 
immediate purpose here, the important consequence of this is that the Gram-Charlier series is not a 
proper asymptotic expansion for piv). In an asymptotic expansion, the remainder term should be of 
higher order than the last term retained. However, if we truncated the series (12) at the (^("^^ term, 
which is 0((T^), we would miss another 0{ct^) contribution coming from cq (equation (15)). In order 
to deal with this problem, let us rearrange the Gram-Charlier expansion by collecting all terms with 
the same powers of a. The result is the so-called Edgeworth series, with the first few terms given by 



p(z/) = (^(z/) 



Z/ (7 



(16) 



Cramer (1946) lists the Edgeworth series to higher order, and he proves that it is a proper asymptotic 
expansion. This proof is directly relevant to our purposes, since it implies that there are no additional 
O(a^) terms hiding in the Gram-Charlier series at m > 6. 

Now we can see the attractiveness of the Edgeworth series for describing the gravitational evolution 
of gaussian fiuctuations: it becomes a series expansion for the evolving PDF in powers of the r.m.s. 
fiuctuation a. This makes physical sense because the Edgeworth series provides an expansion about 
a gaussian probability distribution. If the initial fiuctuations are gaussian, then we expect the terms 
describing successively larger departures from a gaussian PDF to come in with successively higher 
powers of a. For similar reasons, the Edgeworth expansion has recently found applications in stellar 
dynamics as a description of galaxy line profiles (e.g. van de Van der Marel & Franx 1993; Gerhard 
1993). A multivariate Edgeworth expansion is also well-known in kinetic theory: this is the so-called 
Grad solution of the Boltzmann equation (Klimontovich 1982). 

Given equation (16), we can compute the Edgeworth approximation to the PDF provided that we 
can compute Sm to the required order. In this paper we will make use of the second-order approxi- 
mation, 

— SsaHsii') 



1 



3! 



and the third-order approximation. 



p{v) 



1 



j^S^aH^iv) + ^^S^a^H^iv) 



10 



6! 



cp{v). 



(17) 



:i8) 



Although equation (17) contains only a single explicit power of a, it is appropriately described as a 
second-order approximation because the parameter 6*3 remains zero until second order in perturbation 
theory. Similarly, equation (18) is a third-order approximation because 6*4 remains zero until third 
order. 

Of course, if the variable we are interested in is 9 instead of 6, one just adopts the appropriate 
analogues of of Sm-, ^-c, 6*3 — > T3 = (6'^) {O"^) , and so on. Similarly, if we are interested in the 
PDF shape in redshift space, one just needs to use the corresponding analogue of 5^, 5*^5 ^^e §2.7.2). 
The Edgeworth series may also be used to relate 6*3 to other measures of asymmetry like (^|^|) which, 
according to Nusser feDekel (1993), may offer better signal to noise ratio than {6'^) when applied to 
real galaxy surveys (by being less sensitive to the tail of the PDF, i.e., to rare events). This latter 
quantity is not easy to compute directly by perturbation theory (see the appendix of [28]), but it is 
trivial to obtain 



m) = \l-^S^a^ + 0{a'). 
97r 



(19) 



by using the approximation (18) 



2.6 Validity of Perturbation Theory 



The previous results assume that the system never gets to be strongly non-linear, i.e., ci^ ^ 1 is true 
at every scale. But in practice we know that the observed density contrast is very large at small scales. 
One only gets ai <^ 1 for i ^ 8h~^ Mpc. It is thus by no means obvious that there is any range today 
for which perturbation theory might be applicable. On the other hand, it has long ago been noticed 
that linear perturbation theory yields a good description of the variance [or the 2-body correlation 
function ^(t) oc D(ty / D(ti) ^(ti)] measured in fully non-linear numerical simulation. Actually, the 
agreement remains true even for values of the variance as large as ~ 2! 
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Figure 1: The measured S'iii) in a CDM simulation (triangles) is compared to the theory (circles; 
the error bars come form the numerical uncertainty in evaluating the integrals giving 6*3). The left 
panel corresponds to a spherical top-hat smoothing, while the right one corresponds to a gaussian 
smoothing. The vertical dashes mark the limit between the strongly non-linear regime {a > 1) and 
the weakly non-linear one (cr < 1. Courtesy Colombi 1993 & Juszkiewicz et al. 1993a (Ref. [12] & 
[27]). 

It turns out that what holds true for linear perturbation theory results is also true at higher orders. 
Indeed, our top-hat results for 6*3 have been checked against the simulation results of Efstathiou et 
al. 1988, Bouchet and Hernquist 1992, and Lahav et al. 1993 (Refs. [15], [9], & [31]) in the case of 
scale-free initial spectra. In the case of a gaussian smoothing, we checked the theory with the scale 
free- simulations of Weinberg and Cole 1991. In the CDM case, comparison were made with the results 
in Bouchet and Hernquist 1992 (top-hat case, Ref. [9]), as well as those we obtained for that purpose 
which are displayed in figure 1. In all those cases, the agreement between the perturbation theory and 
N-body experiments was excellent (see [27] for further details). 

An analysis of P(S,i) and P(9,i) can be found in Juszkiewicz et al. 1993b (Ref. [28] from which 
we extracted figure 2). It shows the results for a scale-free, ra = — 1, case analyzed with a gaussian 
smoothing. On the left is shown the evolution of the asymmetry measure M3 = (^■^) and M3 = (^|^|). 
Logarithms are base-10. In the lower panel, the solid line shows the prediction of second-order per- 
turbation theory, M3 = Ssa"^, using the value S3 = 3.47 (appropriate to an ra = —1 power spectrum 
and gaussian smoothing filter). Points show measurements from the density fields of the N-body sim- 
ulations, with smoothing lengths of 2, 4, and 8 cells (circles, triangles, and squares, respectively). For 
closer inspection, the upper panel plots the ratios M^/a'^ (top points) and M^/a^ (bottom points), 
with horizontal lines representing the analytic predictions. Error bars mark the la theoretical uncer- 
tainty, i.e., the run-to-run dispersion of eight independent simulations divided by 7^^^. The right side 
of the figure shows the evolution of —M^g = — {9'^) and —M^g = — {0\9\). To plot the lines, we use the 
appropriate value T3 = —2.19, so there are no free parameters to either of these "fits". It demonstrates 
(see also figure 7) that perturbative and N-body results agree to within the la uncertainty when a is 
small, as expected. 
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Figure 2: See legend in text. Courtesy Juszkiewicz et al. 1993b (Ref. [28]) 

Note also that the N-body results for M3 remain remarkably close to the perturbation theory 
prediction even when ci = 1. A similar applicability range was found in the other studies mentioned 
so far. Although these numerical comparisons are not general analytical proofs, they do show the 
usefulness of perturbation theory for the cases of interest in the cosmological context. 

One can also check when the Edgeworth series and perturbation theory are appropriate to describe 
the overall shape of the PDF, by comparing its predictions with the results of N-body experiments. 
The solid lines of figure 3 show the PDF of the smoothed density fields of matter density (left column) 
and velocity divergence (right column), at three different stages in the simulations used above (to 
produce figure 2). The corresponding rms amplitude {a and ag) are shown in each panel, and the 
variables on the horizontal axes are v = bja and v = 9/(Jg. In all the panels, dotted and dashed lines 
show the approximations corresponding to equations (17) and (18) respectively. The approximations 
work well for Ssa < 1 and |z/(t| < 1 (and equivalent requirements for 9). They begin to break down 
outside of that range, as expected. 

2.7 Lagrangian versus Eulerian Approach 

As was shown above, a lot can be accomplished by using the standard Eulerian perturbative approach 
to gravitational instability. However, perturbative equations of motion are often easier to integrate 
when expressed in Lagrangian coordinates. In our experience, this happened at second order for 
fl ^ 1 (see §2.7.1 below). Another example is the redshift space distortion for (S'^), discussed in §2.7.2. 
Moreover, as (^^) grows with time, at any fixed order in perturbation theory, the Lagrangian approach 
is likely to remain valid longer than the Eulerian approach. This is so because the requirement of small 
Lagrangian displacements and gradients is weaker than the requirement (^^) ^ 1. This idea, which 
motivated the Zel'dovich (1970) approximation (which is simply the first order Lagrangian solution), 
remains valid at higher orders. 

In the following, we just outline the derivation of the perturbative solutions, and focus on the 
redshift distortion problem. We also give examples of comparison between Eulerian and Lagrangian 
theory when (^^) ~ 1. 



2.7.1 Lagrangian Formalism In the Lagrangian approach, the primary object of the analysis is 
the particle trajectory instead of the density contrast. A fluid element (a "particle") is indexed by its 
unperturbed Lagrangian coordinate q and its comoving Eulerian position at time t, x(i), is connected 
to q by a displacement field 

x = q + *(i,q). (20) 

An equation for K = V ■ 'if may be obtained by taking the divergence of the motion equation, which 
involves (once Poisson equation is used) the density contrast. The latter is expressed by requiring that 
mass be conserved, i.e., 

p{yi)-J d^q = p{<\)d^q. 

Thus 6 = — 1, if / stands for the determinant of the Jacobian of the transformation from x and 
q. Of course is then only determined up to the addition of any divergence-free displacement field^. 

As in the Eulerian case, perturbative solutions are obtained by means of an iterative procedure, 
as seems to have been done for the first time by Moutarde et al. (1991; it was further developed by 
Bouchet et al. 1992; 1993, refs.[8] & [11]; it also initiated a parallel effort by Buchert and collaborators, 
see the contribution of Weiss and Buchert in this volume, and references therein). Let us denote by 
and M(™), the m — th order part of the (invariant) scalars 

K =V-* = E/*/,/ 

L =\ E/^m(*/,/*m,m - */,m*m,/) (21) 

M = \det[^,J\ 

where $/,m denotes S'^i/Sqm, i.e., the partial derivative of the / — th component of the displacement 
field with respect to the m — th component of the lagrangian coordinate. Focusing on the fastest 
growing mode solutions, we then find 

K{t,q) = g^it) K(^\t„q) + g2{t) L(^\t„q) 

+ g3a(t)M(^Ht,,q) + gsb(t)L(^Ht,,q) ^ > 

where gi{t) is nothing else than the standard linear Eulerian solution D(t)/ D(ti), and g2, gsa, and ^35 
behave near i7 = 1 as 

r, ~ "^0-2/63^2 ^ ^ "'^0-4/77^3 ^ ^ "'^^ 0-2/35 „3 
52 §1, 6'3a~--" §1, 536 §1- 

Further details, and expressions for fiat models with A 7^ are given in [11]. 

The first (linear) term of (22) is nothing else than Zel'dovich solution, since equation (22) implies 
for a potential movement, 

*(i)(r,q) = 5i(0*«(q), (23) 

if we denote by ^('^(q) the initial displacement field. Also, note that ^^^'^(t) = —V • ^f^^) oc gi(t), i.e., 
the Eulerian linear behavior is recovered. The second order growth rate, g2, yields the expression for 
k(^1) introduced earlier [eq. (3)]. And at this order, the solution is again separable as a growth rate 
times a purely spatial part (denoted by tilde) which is set at the initial time: 

*(2)=52(0*(^)(q), with V.*(2)(q) = i(2)(^.0 = ^E(*!?*e™-*!:IO (24) 

(for a potential movement). 

^When needed, one may lift this indeterminacy by restricting one's attention to potential movements, which must 
satisfy Vx x x = 0. This case is relevant, since vortical perturbations decay with time in linear perturbation theory, a 
consequence of the conservation of angular momentum in an expanding universe. Thus one might consider that such 
solutions will apply, if vorticity is initially present, at later times when it has decayed away. 



2.7.2 Real Space- Redshift Space Mapping In redshift space, the appearance of structures is dis- 
torted by peculiar velocities. At "small" scales, this leads to the "finger of god" effect: the clusters 
are elongated along the line-of-sight due to their internal velocity dispersion. This is an intrinsically 
non-linear effect, and we shall not be concerned with it. At "large" scales, the effect is reversed: the 
coherent infiow leads to a density contrast increase parallel to the line-of sight. Indeed, foreground 
galaxies appear further than they are, while those in the back look closer, both being apparently closer 
to the accreting structure (Sargent & Turner 1977; LSS, §76; Kaiser 1987). 

Let us now calculate (^^^ ^^3^ using the Lagrangian approach. Since the unsmoothed density 
contrast is given by ^ = — I, the first terms of its expansion will be given by ^^^^ = —J^^\ 
^(2) = /(I) 2 _ /(2)_ assume an initially gaussian density field. Thus we shall require that the three 
components of the initial displacement field ^^'^ be independent and gaussian, which will insure that 
the density contrast is also gaussian since it is related by a linear operator to the displacement field 
(S, = -V • *(^)). In that case, variance and skewness are given by 

where all averages on the displacement field are taken with respect to the Lagrangian unperturbed 
coordinate2 q. The previous formulae take into account the fact that (^J^^^"^ = ('^^^^) ~ ('^^"^^) ~ ^' 
which also insures that {6) = 0), as well as <^/(i) = (indeed {^$(1)2™+!^ = Q). 

If we call al the variance of any component i of the gradient field = {^^iV^^ (= gl (e^)^ /3), 

which is also gaussian, we have (^J^^^ = 27af = 3 (^J^^^ . The other term in {S'^) involves the 
product /(I) 2/(2) -^Yhich can readily be estimated since, after the second order solution, we have 
/(2) = (1 + g2/gl)L^^y It follows by development that <^/(i) V^^)^ = 6(1 + 52/5D^i- We thus have 
the remarkably simple result 

Ss = 4-2g2/g! + 0{s% (25) 

which is identical to (6). The first term corresponds to the pure Zel'dovich approximation and had 
been found by Grinstein and Wise (1987). 

Let us now consider the case of spherical coordinates, when distances to the observer would be 
estimated by means of redshift measurements. And let us now denote redshift space measurements 
by the superscript z. The redshift space comoving position of a particle located in r(q) = ax(q) is 

= r/(aH) (with H = d/a, where a is the scale factor, while the dot represents a time derivative). 
The real space perturbative expansion (20) is then replaced by 

= q + [1 + fiit)]grm^'Hq) + [1 + Mt)] g2m^'Hq) + 0{e% 

where we have explicitely used the separability of ^'(^) = ^i(i)^(^)(q) and ^^(2) — ^2(i)^(^)(q) 
(eqs. (23) & (24)). For i7 close to I, the logarithmic derivatives of the growth rates are well approxi- 
mated by 

h = (alg^)dg^lda^n^l\ and h = {a/ g^) dg2l da ^ 2^^'^ (26) 

with /i taken from Peebles (1976), and the second from [II]. In the limit of an infinitely remote 
observer, say along the rs-axis, the observed density constrast 6^ in comoving coordinates is simply 
S^(xi, X2, x^), which amounts to approximate spherical coordinates by cartesian ones. All we have 
to do, then, is to replace everywhere in the calculation of S3 the quantity $3™) = g^ $|j™) by (1 -|- 
/m) 9m ^3™^ (for m = 1 and 2). It follows that"^ 

SI = 6 - + + + + + 1 > f > 0. (27 , 

[2 + (l + /i)f - - 



^ Indeed, = (Q[^{<l}]J)q where the subscript x and q refer (respectively) to averages in Eulerian or Lagrangian 

space. This follows if one assumes, as usual, that ensemble averages (...) are equivalent to averages over space, provided 
the volume is large enough that it can be considered a "fair" sample. 

■^All-symmetry breaking terms, which cannot be evaluated without specifying the initial power spectrum, are gathered 
in £. This result was obtained by using Fourier analysis. 



Further details of the derivation may be found in [11]. Of course, we recover the real space result (25) 
if we set fi = f2 = 0. On the other hand, if i7 = 1, we have = 1 = /2/2 (and §2/ 9\ = —3/7), which 
yields 6*3 = (35 + S)/! while, for i7 = 0.1, 6*3 ~ (34.5 + QAS)/! . Since we used an "infinitely remote 
observer" approximation, the formula (27) strictly applies only in the limit of large volumes. In any 
case, at least in that limit, it clearly shows that the ratio 6*3 is nearly independent of the value of i7, 
nor is it affected by redshift space distortions. 
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Figure 4: See legend in text. Courtesy Hivon et al. 1993 (Ref. [23]) 



We have recently extended those calculations by including smoothing and by dropping the simpli- 
fying assumption of an infinitely remote observer {i.e., without approximating spherical coordinates 
by cartesian ones, but keeping the "large sample limit", as in [29]). The results are displayed in 
figure 4, which is extracted from [23]. They confirm that the skewness of an initially gaussian field is 
essentially the variance squared times a coefficient that depends mostly on smoothing and the initial 
variance. Most of the redshift or i7 dependence is contained in the variance dependence. 

Finally, let us mention that T^^ can be obtained by similar techniques. Indeed, with the same 
notations than above 

Since 6'(^) = figiK^'^\ and 6'(^) = f2g2K^^^ - /i(Vg*(^) • Vg*(^)), the assumption of gaussian initial 
conditions in Fourier space (see note 2) yields 

r3 = 4(i-^%).-^i)--. (28) 

Once the expressions (26) are inserted in this formula, one recovers equation (11) when no smoothing 
correction is included'*. By using Zel'dovich approximation, one would not include the second order 
term, which then yields ~ — 2i7~°'^. When a top hat smoothing is included, Zel'dovich approx- 
imation yields ~ — i7~°'^(ra -|- 1) instead of the correct answer Ts — —fl~^'^{n — 5/7). Even the 
sign of is wrong for — l<ra<5/7!Itis also interesting to note that on scales where one might 
want to measure T3, the effective index of the power spectrum is close to n = —1. Then Zel'dovich 
approximation leads to an essentially unskewed PDF, whose shape is mainly governed by its kurtosis 
(which is also not computed exactly appropriately). 

2.7.3 Approximation of Non-Linear Dynamics So far, we have used Lagrangian perturbative 
solutions at the appropriate order to obtain "exact" results, in the regime when such an approach 
should be applicable. But one may also think of using these solutions as approximation to the real 

*More precisely, T3 ~ — 2Q~'^^^(1 -|- If instead of using eq. (26) one uses limited expansions around 
n = 1, !,e,, /i ~ Q"/^ and /2 ~ 2Q^/^then one obtains a slightly different expression T3 ~ -2Q-^/^(l + f Q-^^^^^)), as 
in [6]. 



dynamics in the not-so-weakly non-linear regime. Indeed, Zel'dovich [43] idea was to use his ballistic 
approximation (where gravitational acceleration is ignored), even when the density contrast S becomes 
large. It lead to the development of pancake theory. Of course, as stated by Zel'dovich, "the analytic 
evaluation of the error is extremely difficult". This lead to many comparisons of Zel'dovich approxi- 
mation with the exact dynamics, either with simulations (e.g., Doroshkevich et al. 1980) or rigorous 
Eulerian perturbative theory (Grinstein and Wise 1986, Juszkiewicz et al. 1993). 
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Figure 5: Computed amplitude 6 (left panels) and velocity field (right panels), for a spherically 
symmetric overdensity (top panels) and an underdensity (bottom panels). The solid line shows the 
exact result. The long dashes correspond to the second-order Lagrangian approximation, the short 
dashes to the first order one, i.e., Zel'dovich approximation, the short dashes-dots to the second-order 
Eulerian approximation, and the dots to Eulerian linear theory. This is for the case i7 = O.I, but the 
results do not depend much on this value. Courtesy Bouchet et al. 1993 (Ref. [11]). 

Despite its limitations, this approximation turns out to be amazingly good, at least when the initial 
field is smooth enough. Indeed, it is widely used today^, to the point of being employed to address 

^For instance to predict the weakly non-linear evolution of the moments of the one point probability distribution 
function of the density field (Betancort-Rijo 1991; see also Hoffman 1987 for the variance only), or of the distribution 
itself (Kofman et al. 1993, Padmanabhan 1993), or in reconstruction methods to recover the "initial conditions" from 
present day observations {e.g., Nusser et al. 1991, Nusser and Dekel 1992, Gramman 1992, Lachieze-Rey 1993a). 



questions where it is not really appropriate. The success of Zel'dovich approximation brought about 
many attempts to do better by correcting its shortcomings, and a number of them were discussed in 
this meeting. The question of how good are those approximations to describe the exact dynamics thus 
naturally arises. Here we restrict ourselves to rigorous perturbative expansion. 

We have compared in [11] the first and second order solutions in Eulerian and Lagrangian pertur- 
bation theory to spherically symmetric cases whose evolution is analytically known. For instance, the 
results of those approximations were checked for the density and the divergence of the velocity field 
in the spherical top-hat case, when its amplitude is varied. The figure 5 shows the result of another 
comparison, when the profile is smooth. These lead us to the following approximate ranking (at least 
for moderate final density contrasts ~ 1): 

density contrast : 

Lagrangian second order > Zeldovich ^ 
Eulerian second order > Eulerian linear theory. 

Velocity field : 

Lagrangian second order ~ Eulerian second order > 
Zeldovich ~ Eulerian linear theory. 

Here the signs ">" and "~" mean respectively "more accurate than" and "of comparable accuracy 
to". For relatively large final density contrasts, the Eulerian approach becomes particularly inefficient, 
except for the velocity field, for which it tends however to be less accurate than the Lagrangian one. 
The second order Lagrangian approach gives, for moderate final S, an excellent approximation of the 
density contrast and the velocity field. Its seems to be able to reproduce density contrasts as large as 
ten. 

Direct comparisons with numerical simulations (Ref. [11]) appear to confirm that ordering. In 
essence we find that second-order Lagrangian perturbation theory works well for density constrast up 
to unity. We expect only little improvements to be brought at higher orders. 



3 Application to Data 

3.1 ffiAS Density Field 

Bouchet et al. (Refs.[7] & [10]) have measured the count probability distribution function (CPDF) 
in a series of 10 volume-limited sub-samples of a deep redshift survey of IRAS galaxies^. Counts 
were performed in spherical cells, using redshifts as distance estimators. Among other things, they 
deduced S3 by computing various centered moments of the CPDF. Figure 6 shows on the left an 
equal weight average in bins of values of log (^^) (,^2 is just an equivalent notation in the absence of 
discreteness corrections''), the average value S3 = 1.5 ± 0.5 (dashes), as well as the values inferred 
from measurements of Q in the non-linear regime from optical data (triangles), from measurements of 
skewness and variance on the QDOT sample (error bars on left) and the theoretical prediction (solid 
line on left) from perturbation theory for a power spectrum of index n = —1.4 and bias bi = 1 (se 
below). 

But since Strauss et al. (1992a) show that the galaxy densities in cores of clusters determined from 
IRAS galaxies are systematically lower than those determined from optically selected galaxies, they 
also did counts were galaxies associated with cluster cores were assigned an extra weight corresponding 

^The sample consists of 5304 galaxies with 60 micron flux density above 1.2 Jy, selected over 87.6% of the sky. The 
selection criteria for the galaxies are given in Strauss et al. (1990) and Fisher (1992), and the data for the brighter half of 
the sample are given in Strauss et al. (19926). IRAS galaxies are a dilute tracer of the galaxian density field (Strauss et 
al. 1992a), with typically 1/3 the number density of galaxies appearing in optically selected samples of comparable depth. 
Thus one can explore only the low-density limit, but the large volume covered by this sample allows many independent 
volumes of a given size at a given number density to be probed. 

^Otherwise, (S^) = 1/1V-F6- 



to the ratio of the optical and IRAS density estimates. They refer to those counts as the boosted 
counts. As the right panel of figure 6 shows, it makes a quite noticeable difference. In particular, the 
average value of S3 is appreciably higher (for 0.1 < ^2 < 10), namely 3.71 ± 0.95, which is more than 
twice the value found for the unboosted case. This naturally raises the issue of biasing. 
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Figure 6: See legend in text. Courtesy Bouchet et al. 1993 (Ref. [10]) 

Indeed, perturbation theory tells us the S3 to expect from the gravitational evolution of gaussian 
initial conditions. But it is widely believed that the density field traced by galaxies selected in some 
way (in the optical, or the infrared, etc. ), ^g, may not be identical to the S considered so far. Let us 
suppose, as in [28], that the smoothed galaxy field Sg is a local, but non necessarily linear function of 
Si, i.e., Sg = B(Si). Then, by using a Taylor expansion, 

C - (^^') _ ^3 , 

with hi = B'(0) and 62 = B"(0). This shows that such a biasing preserves the proportionality of the 
skewness to the square of the variance®. One can therefore conclude that the data is compatible with 
the hypotheses of gaussian primordial fiuctuations and local biasing. Furthermore, as mentionned 
earlier [eq. (4)], a non-gaussian field would yield a term oc a, the proportionality depending of course 
on the initial value of the skewness (e"^) (it would also introduce a term involving the initial kurtosis, 
in the case of a non-linear biasing). One will thus need to specify a particular non-gaussian model 
(e.g., cosmic strings or textures which are considered in [12] & [3] respectively) to assess whether the 
data puts a strong constraint on such a model. 

3.2 POTENT Velocity Field 

The POTENT Velocity Field may be used to check whether relation (11) holds. Since IRAS suggests 
that the hypothesis of gaussian fiuctuations is compatible with the data, one can attempt to use (11) 
to constraint the value of i7. Figure 7 shows the measured value of Ts in POTENT samples of various 
volumes (but at a fixed scale). Small volumes correspond to the best data, while larger volumes offer 
better statistics. The "best" determination is probably for a volume of radius around 4500 km/s. In 
order to further assess the effect of sampling, observational errors, finite volume sizes, etc. , similar 

*Fry & Gaztanaga (1993) have independently derived equation (29), and they have generalized the result in an 
important way: by expanding 8g = f{8) in higher order Taylor series, they show that a local biasing function preserves 
all of the moment relations predicted by perturbation theory, in the limit of small fluctuation amplitude. 



(29) 



measurements were performed on a series of CDM simulations with 2 different values of i7 (solid 
squares correspond to i7 = f , while solid triangle correspond to i7 = 0.3. 

Taken at face value, the measurement are consistent (within one expected standard deviation) 
with i7 = f , while i7 = 0.3 would be excluded at the level of two standard deviations. It should be 
emphasized, though, that the data analyzed is still preliminary (i.e., it corresponds to today's data, 
but before it was cleaned, calibrated, and put together properly). Furthermore, the real statistical 
significance of these "standard deviations" is not really known - the present data is too small to 
estimate the possible sampling errors due to missing large scale power. One may also worry about the 
effect of imperfectly correcting for an inhomogeneous Malmquist bias. Our estimate off) is thus more 
of an illustration of the fact that the method should be taken seriously and should have a practical 
value when larger datasets become available. Future work (in particular concerning the effect of an 
only partially corrected inhomogeneous Malmquist bias) will tell us how to treat these indications. In 
any case, it provides an example showing that this method can work. 
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Figure 7: Measurements of — Ts are denoted by solid circles; the attached error bars correspond 
to a series of measurements made with "fake" data obtained by offsetting the real measurements 
by some amount given by the expected size of the observational errors, in the standard POTENT 
analysis manner; the open circles correspond to the raw measurements. See legend in text for further 
comments. Courtesy Bernardeau et al. 1993 (Ref. [6]). 



4 Summary and Discussion 

In this paper, we have described the results of our studies of the relationship between statistics and 
dynamics of a weakly non-linear, self-gravitating, pressure-free fiuid in an expanding Universe. In 
particular, we have investigated how gravitational instability drives the distribution of density and 



velocity fluctuations away from the initial state, which we assume gaussian for the most part. 

We compared our analytical, perturbative calculations of low-order reduced moments of PDFs of 
the density and velocity divergence held with N-body simulations. We found excellent agreement over 
a surprisingly wide dynamical range, all the way to ci ~ 1 and (Tg ~ 1, when the perturbative series is 
expected to blow up. Using the Edgeworth expansion, we have calculated the PDFs of S and 9. These 
expansions truncated after the flrst non-trivial terms turn out to give a good approximation to the 
real thing provided that |z/(t| < 1 and Ssa < 1 (and similar requirements for ag). 

We have investigated the effects of redshift space distortion on the lowest order non trivial cumulant 
of the PDF of S - the skewness (S'^) - and found that while both skewness and variance are modifled 
by the distortion, the ratio S3 = {S'^) {S"^) ^ is little affected, at least in the weakly non-linear regime 
(a ^ 1). We have also investigated the possible effects of biasing of the distribution of galaxies with 
respect to that of mass. We found that when the biasing is a local function of S, the scaling (S'^) oc a"^ 
(expected from gravitationally induced skewness of an initially gaussian held) remains valid for counts 
of galaxies. These two effects - the lack of redshift distortion and the preservation of the scaling law - 
can make S3 measurements a powerful tool in distinguishing gaussian initial conditions from strongly 
non gaussian alternatives, when one expects (^■^) oc a'^. 

The sensitivity of (6'^) to fl can be used to measure the density parameter from the peculiar velocity 
data alone rather than from heterogeneous data sets (like in the IRAS-POTENT comparison). Such 
a measurement is also insensitive to whether or not galaxies trace mass. The only requirement is 
that they do trace the true peculiar velocity held. Estimating i7 from (^■^) is somewhat similar to 
the "reconstruction method" of Nusser & Dekel (1993). An important difference between our two 
approaches is that Nusser & Dekel use Zel'dovich approximation while we use rigorous perturbation 
theory. The problem with Zel'dovich approximation is that it fails rather badly in recovering the true 
numerical values of the reduced moments of the PDFs (in particular those of the smoothed velocity 
held and/or 9, as was briefly discussed in §2.7.2). Since we showed that the estimate of i7 depends 
precisely on these moments, one wonders how accurate can be the resulting estimates of i7. There 
is also an observational problem, common to both approaches. The currently available sample may 
not be deep enough yet (as always seems to be the case in cosmology). The coherence length of the 
flow appears comparable to the size of the region mapped. Such a sample thus does not appear to be 
"fair"enough to reliably measure the true distribution of 9, or even its flrst three moments. 

There is clearly more work to be done in the area of perturbation theory, in particular to assess the 
combined effects of weakly non-gaussian initial conditions, a non-linear bias, for density observations 
in redshift space. Similarly, further studies are required to check whether perverse effects (flniteness of 
the sample, inhomogeneous Malmquist bias, etc. ) might spoil the simple method using T3 to measure 
i7. But, even if these problems taken together turn out to be too formidable to be fully analytically 
tractable, they are now posed well enough, and the interest of doing so sufficiently well established, 
that they may be resolved by Monte Carlo analysis. 

The ultimate goal of our efforts, presented here, is to use the PDFs of S and 9 to test for a generic 
family of models, based on gravitational instability acting on gaussian initial conditions. So far, we 
have reliably tested our perturbative calculations with N-body experiments. There are indications of 
agreement with real data as well. Hopefully, in the near future, new generations of catalogs, like the 
Digital Sky Survey, will flnally tell us whether all we are looking at is just a gaussian random process 
modifled by plain gravity. 
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